Multiscale Random- Walk Algorithm for Simulating Interfacial Pattern Formation 
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We present a novel computational method to simulate accurately a wide range of interfacial patterns 
whose growth is limited by a large scale diffusion field. To illustrate the computational power of 
this method, we demonstrate that it can be used to simulate three-dimensional dendritic growth in 
a previously unreachable range of low undercoolings that is of direct experimental relevance. 
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Interfacial patterns form spontaneously in a wide range 
of physical and biological systems where the motion of 
an interface is limited by one (or several) diffusion fields, 
each one obeying the diffusion equation 



dtu = ZJV^u 



(1) 



with specified boundary conditions on the interface. 
Classic examples include various cellular, dendritic, or 
eutectic solidification patterns (where u represents the 
temperature or an impurity concentration) pi, dendrite- 
like branched patterns formed during electrochemical de- 
position (where u is some ion concentration) and com- 
plex growth morphologies produced by bacterial colonies 
under stress (where u can represent the concentration of 
some nutrient or a signaling agent) ^. 

When attempting to accurately simulate the growth of 
such structures, one is generally faced with two major dif- 
ficulties. The first one is front tracking, which requires 
to resolve accurately the boundary conditions imposed 
on u at the evolving interface. Dendritic solidification, 
where even a weak crystalline anisotropy crucially influ- 
ences the morphological development epitomizes this dif- 
ficulty. The second problem is the large disparity of scale 
between the growing structure and the diffusion field sur- 
rounding it. The physical origin of this disparity is es- 
sentially dimensional. The diffusion field decays ahead of 
the growth structure on a length scale I ~ D /v, where 
V is the velocity of the advancing interface, whereas the 
characteristic scale of the structure, e.g. the tip radius p 
of a growing dendrite, is itself the geometric mean of I and 
a short length scale cutoff proportional to the interface 
thickness. Thus, for small growth rate, I can be several 
orders of magnitude larger than p, and simulations that 
resolve simultaneously the details of the interfacial pat- 
tern and the diffusion field become extremely difficult. 

Whereas various methods have been successfully devel- 
oped to handle front tracking, bridging the length scale 
gap between p and I has remained a major computational 
challenge. A natural idea to overcome this problem is 
to use multigrid or adaptive mesh refinement algorithms 
that make the grid progressively coarser away from the 
interface [0-pl. However, such methods need to dynam- 
ically adapt their grids to follow the moving interface. 
This is a nontrivial task, and quantitative simulations 



of dendritic crystal growth at low undercooling have re- 
mained restricted to two dimensions (2-d) g. 

In this letter, we present a novel hybrid computational 
approach that efficiently bridges this length scale gap. 
Over most of the computational domain, the diffusion 
equation is simulated by an ensemble of off-lattice ran- 
dom walkers that take longer, and concomitantly rarer, 
steps with increasing distance away from the growing in- 
terface. This drastically reduces the computational cost 
for evolving the large-scale field. Moreover, a short dis- 
tance away from the interface, this stochastic evolution 
is connected to a finite-difference deterministic solution 
of the interface evolution. This conversion, in turn, re- 
duces the inherent noise of the stochastic method to a 
negligibly small level at the interface. This approach is 
relatively simple to implement in both 2-d and 3-d while 
being at the same time quantitatively accurate. Here we 
sketch the method and then report results that demon- 
strate its capability to yield new quantitative predictions 
testable by experiments in the context of dendritic crystal 
growth. For clarity, we expose the method in this context 
although it will become clear below that it is general. 

Let us consider a solid-liquid interface whose motion is 
limited by heat diffusion and define a scaled temperature 
field u that is zero in equilibrium and equal to —A in the 
liquid far from the interface, where A is the dimensionless 
undercooling. At the interface, u satisfies the well-known 
boundary conditions 
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(2) 
(3) 



i=l 



corresponding to heat conservation and local thermody- 
namic equilibrium at the interface, respectively, where 
Vn is the normal velocity of the interface, do is a micro- 
scopic capillary length, 9i are the local angles between 
the normal fi to the interface and the two local princi- 
pal directions on the interface, Ki are the principal cur- 
vatures, and the function a{h) describes the orientation 
dependence of the surface energy. 

The basic idea of the method is to divide space into 
an 'inner' and an 'outer' domain as illustrated in Fig. n. 
The inner domain consists of the growing structure and a 
thin 'buffer layer' of liquid surrounding the interface. The 



outer domain corresponds to the rest of the hquid, and 
is much larger than the inner domain at low undercool- 
ing. In the inner domain, we solve deterministically the 
diffusion equation on a fine uniform mesh. Moreover, for 
the present crystal growth application, we handle front 
tracking using a phase- field approach 0], and time-step 
explicitly both u and the phase-field in the inner region 
using the same procedure as Karma and Rappel |^ . The 
geometry of Fig. ffl, however, implies that any other front 
tracking method that solves the diffusion equation on a 
uniform mesh could be used instead. Moreover, since the 
boundary conditions on u need not necessarily be those 
defined by Eqs. || and ^, the method obviously extends 
to other diffusion- limited pattern forming systems. 

In the outer domain, the diffusion equation is simulated 
stochastically by an ensemble of off-lattice random walk- 
ers. The idea of solving the diffusion or Laplace equation 
with an ensemble of random walkers is well-known and 
has been used previously to simulate diffusion-limited 
growth [^ and Hele-Shaw flow |l^. The main new fea- 
ture of our method is that we have separated the solid- 
liquid interface from the boundary at which the conver- 
sion from the deterministic to the stochastic solution of 
the diffusion equation takes place. This separation yields 
two essential benefits. Firstly, it makes it possible to 
use the phase-field approach with its proven accuracy 
to simulate the interface evolution and to resolve even 
a weak crystalline anisotropy, without being affected by 
the details of the conversion process. Secondly, in the 
buffer layer between the solid-liquid interface and the 
conversion boundary, the temperature obeys the deter- 
ministic diffusion equation. Consequently, the noise cre- 
ated by the stochastic release and impingement of walkers 
is rapidly damped away from the conversion boundary. 
Hence, the amplitude of temperature fiuctuations at the 
solid-liquid interface can be reduced to an insignificant 
level without much cost in computation time by increas- 
ing the thickness of the buffer layer. 

To connect the inner (deterministic) and outer 
(stochastic) solutions, we have to supply a boundary con- 
dition for the integration of the inner region, and we 
must specify how walkers are created and absorbed at 
the boundary between inner and outer regions. Both pro- 
cesses are handled by using a coarse-grained grid that 
is superimposed on the fine grid of the inner region as 
shown in Fig. H, Cells of the coarse grid on the border 
between inner and outer region, shaded in Fig. n, are 
called conversion cells. The temperature in a conversion 
cell is related to the local density of walkers, 



-A(l-m,(t)/M) 



(4) 



where rriiit) is the number of walkers in conversion cell 
number i at time t, and M 3> 1 is a fixed integer. 
Hence, an empty cell corresponds to u = —A (the initial 
state), whereas a box containing M walkers corresponds 
to w = 0. This formula is used in each timestep to ob- 
tain the boundary condition for the integration of the 




. • • • 

• • • 

.. • 

1 • • . . • 

— — — - I I I I I I I I 1 

I • 
* • • 

"^ • • 

'__^ — ^ 



FIG. 1. Top: snapshot of a 2-d simulation for A = 0.1 and 
£4 — 0.025 showing the solid-liquid interface (solid line), the 
inner-outer domain conversion boundary (dashed-line) , and 
the random walkers (dots). Only a small part of the outer 
domain and one out of 50 walkers are shown for clarity. Bot- 
tom: enlarged view of the conversion boundary (dashed line) 
showing the fine and coarse grids. Shaded cells are conversion 
cells and walkers (dots) are restricted to the outer region. 



inner region. Next, we determine the quantity of heat 
that fiows in or out of each conversion cell from the in- 
ner region, and add this amount to a reservoir variable 
Hi (t) which describes the heat content of conversion cell 
number i. If this variable exceeds a critical value He, a 
walker is created and H^ is subtracted from the reservoir. 
Conversely, if Hi{t) falls below —Hc,b, walker is absorbed 
and He is added to the reservoir. This procedure assures 
that walkers are created and absorbed at a rate which 
is proportional to the local heat fiux, and each walker 
corresponds to the same discrete amount of heat. 

Evidently, as the structure grows the geometry of the 
conversion boundary and the configuration of the coarse 
grid need to be periodically updated in order to maintain 
a constant thickness of the liquid buffer layer. This pro- 
cedure, however, is straightforward since the structure of 
the grids does not change. 

In the outer region, each walker is represented by a set 
of variables indicating its position and the time it has 
next to be updated. To update a walker, a new position 
is randomly selected with a probability distribution given 
by the diffusion kernel. 



P{x',t'\x,t) 
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4L>(i' - t) ' 



(5) 



where x and af' are the old and new positions of the 
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FIG. 2. Result of a 3-d dendrite growth simulation for 
A = 0.05 and 64 = 0.025 with snapshots of the 3-d structure 
(top) at the times corresponding to the arrows. This run took 
6 hours on 64 processors of the CRAY T3E at NERSC and 
used up to 5 X 10^ walkers. 



walker, respectively, and t' ~ t is the time increment 
between updates. This representation of the diffusion 
equation is widely used in quantum Monte Carlo meth- 
ods ||ll|] . The key improvement that makes the algorithm 
efficient in the present context is the introduction of a 
variable step size: we allow walkers to take progressively 
larger steps with increasing distance away from the in- 
terface, and to be concomitantly updated more rarely, 
which does not affect the quality of the solution near the 
solid-liquid interface. Adaptive steps have been previ- 
ously used to speed up simulations of diffusion-limited 
aggregation, albeit in a simpler Laplacian context where 
the stepping time is irrelevant and only one walker at 
a time is simulated P4 ]. Typically, we vary the aver- 
age step size between a value comparable to the spacing 
of the inner mesh to about 100 times that value. Our 
test dendrite computations show that the far field can be 
evolved at essentially no extra cost: the program spends 
most of its time in the inner region and for the walkers 
which are near to the conversion boundary and have to 
make small steps. Thus this 'adaptive step' implemen- 
tation yields essentially the same benefits as an adaptive 
meshing algorithm, while avoiding the overhead of re- 
gridding. Finally, our method can be easily parallelized 
as will be discussed in more details elsewhere. 

To illustrate our method, we focus here on the initial 
stage of dendritic solidification, during which four (six) 
primary arms in 2-d (3-d) emerge from a structureless nu- 
cleus, as shown in the example 3-d run of Fig. 0. While 



this transient regime has recently been investigated nu- 
merically in 2-d and experimentally in 3-d |13| , it has not 
yet been explored by simulations in 3-d. 

The present simulations started from a small spheri- 
cal solid nucleus with a uniformly undercooled tempera- 
ture M = — A, and fully exploited a cubic symmetry (i.e. 
a{n) = {l-3ei)[l-4:ei/{l-3ei){ni + n^ + ni)] with the 
Cartesian axes defined parallel to the [100] directions) to 
reduce simulation time. The value e^ = 0.025 that corre- 
sponds to the experimentally estimated anisotropy value 
for pivalic acid (PVA) jl4| was used in all the simula- 
tions reported here. Results for other anisotropics and 
that pertain to the 3-d morphology of the dendrite tip 
will be discussed elsewhere. To obtain quantitative data 
on the growth transients, we recorded the arm length 
L(i), the tip velocity v{t) = L{t), the tip radius of cur- 
vature p{t), and the total volume of solid Vs{t). In Fig. 
0, we show v{t) and p{t) for a 3-d run at A = 0.05 to- 
gether with the time-dependent tip selection parameter 
a*{t) — 2df)D/[p'^{t)v{t)] and the steady-state velocity 
Vss calculated by a boundary integral method. Two re- 
sults are particularly noteworthy. Firstly, (T*{t) becomes 
essentially constant as soon as the arms have emerged 
from the spherical seed, whereas both v(t) and p{t) are 
far from their steady-state values. Physically, this is a 
direct consequence of the fact that v{t) and p{t) evolve 
slowly on the tip diffusion time scale p^/D where a* is 
established. Secondly, we find that, as the undercooling 
is lowered, the volume of the dendrite (or its area in 2- 
d) approaches that of a sphere (circle) growing at the 
same undercooling. The latter is readily obtained from 
Zener's well-known similarity solution, which yields that 
the radius of a d-dimensional sphere grows as ^/ApdDt, 
where the Peclet number pd(A) is implicitly defined by 
A = p;j/%xp(pd) /~ s-'^/^e-^ds. Both in 2-d and 3-d, 
the volume of the dendrite grows slightly faster than the 
one of the sphere, but for the lowest undercoolings we 
could attain, the final volume differed by only 20% from 
this prediction, even though the arms were already very 
well developed. 

The above observations are in good agreement with 
theoretical expectations for 2-d growth transients. In 
particular, Almgren et al. have analyzed the related 
problem of anisotropic Hele-Shaw flow (i.e. Laplacian 
growth) at constant flux |l5[ . Using an exact solution for 
the Laplacian field around a cross and exploiting the con- 
stancy of a* , they constructed a self-afhne scaling shape 
for the arms. The length and width of this shape grows 
as t"' and t^ , respectively, with a — 3/5 and (i = 2/5. As 
subsequently remarked by Brener ||l6|| , the diffusion equa- 
tion can be replaced by Laplace's equation on the scale 
of the dendrite as long as L/^/dI <C 1, and hence 2-d 
growth transients should obey this scaling with a flux set 
by the diffusive far field of Zener's similarity solution in 
the low undercooling limit. In Fig. ya, we plot the func- 
tions a{t) ^ d{\nL)/d{\nt) and v{t) = d{\nVs) / d{\nt) . 
For an exact Laplacian scaling in 2-d, a{t) = 3/5 and 



i/(t) = 1. With decreasing undercooling both curves be- 
come flatter and indeed approach the expected Laplacian 
scaling. The slow rise with time of both curves can be 
attributed to diffusive corrections to the Laplacian scal- 
ing due to the slow increase in time of L/\/lji. Recently, 
Provatas et al. have reported scaling exponents that dif- 
fer from the Laplacian prediction and that appear to be 
independent of A for small A [Q. We note, however, 
that they used the distance from the tip to the time- 
dependent base (where the dendrite shaft is narrowest) 
to scale their results instead of L{t). Since L{t) is the 
only relevant scaling length for both the shape and the 
diffusion field in the Laplacian limit, we believe that these 
exponents are spurious. When L{t) is used as a scaling 
length, our results are consistent with an approach to 
Laplacian scaling in the limit of vanishing undercooling. 

It is simple to heuristically generalize some of the above 
ideas to 3-d. If we assume, as a reasonable first approx- 
imation, that the arm shape is axisymmetric and has a 
scaling form r{x,t) = t^f (x/f^), where x is the growth 
direction, the constancy of cr* imposes that 4/3— a— 1 = 0. 
Furthermore, assuming that the volume of the dendrite 
grows approximately as the one of the 3-d similarity so- 
lution, we have in addition v — a + 2f] — 3/2. These 
two conditions yield a — 2/3 and f3 — 5/12. Fig. ||b 
shows the functions a{t) and iy{t), defined as before, in 
3-d. As in 2-d, the curves approach the predicted expo- 
nents with decreasing undercooling, but the differences 
remain larger than in 2-d even for the lowest undercool- 
ing. This can be partly accounted for by the fact that the 
tip velocity, and hence also L/i/Dt, is much larger in 3-d 
than in 2-d at equal undercoolings. We must emphasize 
that in the absence of an exact 3-d Laplacian solution, 
and in view of the above assumptions, no claim is made 
here that these 3-d exponents are exact or that a scaling 
regime exists asymptotically at small undercooling in 3- 
d. We content ourselves with the fact that they describe 
reasonably well our present simulations. 

In conclusion, we have presented a novel computational 
approach that can resolve accurately the details of a com- 
plex branched structure and its large scale surrounding 
diffusion field. The method can be combined with many 
of the existing front tracking methods and should be ap- 
plicable to a wide range of diffusion- limited pattern form- 
ing systems. Furthermore, we have demonstrated its fea- 
sibility in the non-trivial test case of dendritic growth in 
a range of parameters previously unreachable in 3-d. 
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FIG. 3. Functions a{t) and i/(t) vs scaled time t/to{A) in 
2-d (a) and 3-d (b). In order to show all curves on the same 
plot, to (A) was defined to be the time at which the tips are 
ahead of the grooves by about two tip radii. In order of de- 
creasing A, toD/dl = 2.29 X 10*, 9.03 x 10^ 4.36 x 10^ 
5.35 X 10^ and 1.39 x 10^" in 2-d, and 4.18 x 10^ 5.58 x 10^ 
and, 1.17 x lO" in 3-d. 
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